Load in libraries

library(tidyverse)
library(readxl)
library(mgcv)
library(emmeans)
library(ggeffects)
library(lme4)
library(formattable)
library(ggpubr)
library(pzfx)
library(DT)
library(cowplot)
source("scripts/utils/utils.R")
modernacolors()
## <environment: R_GlobalEnv>

Load data

# load data
d <- read_excel("processed_data/ELISA.xlsx")

Observed data visualization

Supplementary Figure 3. Impact of mRNA-1273 Prime-boost Interval on S2-P-specific Serum Binding IgG Antibody Titers Through 24 Weeks Post-boost

#-------------------------------------------------------------------------------
# Supplementary Figure 1
#-------------------------------------------------------------------------------

# color panel
panel0 <-
  c(
    modernadarkblue2,
    modernaorange,
    modernapurple,
    modernagreen,
    modernablue,
    modernared,
    modernadarkgray
  )
dosepanel <- c("#cb6939", "#8264cb")
LLOQ <- 25  # lower limit of quantificatin


# trend over time
raw1 <-
  d %>%
  filter(interval != "PBS") %>%
  mutate(interval = factor(interval)) %>%
  mutate(
    .,
    interval = forcats::fct_recode(
      interval,
      "8-week interval" = "8",
      "6-week interval" = "6",
      "4-week interval" = "4",
      "3-week interval" = "3",
      "2-week interval" = "2",
      "1-week interval" = "1",
      "Prime only" = "0"
    )
  ) %>%
  mutate(., interval = factor(.$interval, levels = levels(.$interval)[8:1])) %>%
  mutate(
    day = forcats::fct_recode(
      factor(day),
      "Prime only" = "56",
      "1 wk post-boost" = "64",
      "2 wk post-boost" = "73",
      "4 wk post-boost" = "87",
      "8 wk post-boost" = "115",
      "12 wk post-boost" = "143",
      "16 wk post-boost" = "171",
      "20 wk post-boost" = "199",
      "24 wk post-boost" = "226"
    )
  ) %>%
  mutate(dose = forcats::fct_recode(
    factor(dose),
    "mRNA-1273 1 µg" = "1",
    "mRNA-1273 10 µg" = "10"
  )) %>%
  ggplot(aes(
    x = factor(day),
    y = value,
    color = factor(interval)
  )) +
  geom_boxplot(outlier.shape = NA) +
  geom_jitter(position = position_jitterdodge(jitter.width = .2)) +
  geom_hline(aes(yintercept = log10(LLOQ), lty = "LLOQ")) +
  scale_linetype_manual(values = 2) +
  facet_wrap( ~ dose, nrow = 3) +
  theme_bw() +  theme(panel.grid.major = element_line(colour = "grey100")) +
  labs(
    x = " ",
    y = expression(log[10] ~ S2P ~ IgG ~ titer),
    color = "",
    linetype = ""
  ) +
  scale_color_manual(
    values =   c(
      modernared,
      modernablue,
      modernagreen,
      modernapurple,
      modernaorange,
      modernadarkblue2,
      modernadarkgray
    )
  ) +
  scale_y_continuous(breaks = 1:7) +
  theme(
    axis.text = element_text(size = 12),
    axis.title = element_text(size = 16),
    strip.text = element_text(size = 14),
    legend.text = element_text(size = 14),
    legend.title = element_text(size = 16)
  )

# animal-level line plot
raw2 <-
  d %>%
  filter(interval != "PBS") %>%
  mutate(
    interval = factor(interval),
    animal_id = paste("animal", animal_id),
    animal_id = factor(
      animal_id,
      levels = c(
        "animal 1",
        "animal 2",
        "animal 3",
        "animal 4",
        "animal 5",
        "animal 6",
        "animal 7",
        "animal 8",
        "animal 9",
        "animal 10"
      )
    )
  ) %>%
  mutate(
    .,
    interval = forcats::fct_recode(
      interval,
      "8-week interval" = "8",
      "6-week interval" = "6",
      "4-week interval" = "4",
      "3-week interval" = "3",
      "2-week interval" = "2",
      "1-week interval" = "1",
      "Prime only" = "0"
    )
  ) %>%
  mutate(., interval = factor(.$interval, levels = levels(.$interval)[8:1])) %>%
  mutate(dose = forcats::fct_recode(
    factor(dose),
    "mRNA-1273 1 µg" = "1",
    "mRNA-1273 10 µg" = "10"
  )) %>%
  ggplot(aes(
    x = day,
    y = value,
    color = factor(animal_id)
  )) +
  geom_line() +
  geom_jitter(position = position_jitterdodge(jitter.width = .2)) +
  geom_hline(aes(yintercept = log10(LLOQ), lty = "LLOQ")) +
  scale_linetype_manual(values = 2) +
  facet_grid(dose ~ interval) +
  theme_bw() +  theme(panel.grid.major = element_line(colour = "grey100")) +
  labs(
    x = " ",
    y = expression(log[10] ~ S2P ~ IgG ~ titer),
    color = "",
    linetype = ""
  ) +
  scale_color_manual(
    values = c(
      "#c25dba",
      "#72b24a",
      "#7561cf",
      "#ce9b44",
      "#807dc3",
      "#7a7d37",
      "#4faad9",
      "#c95c3f",
      "#4bac88",
      "#c8577b"
    )
  ) +
  scale_y_continuous(breaks = 1:7) +
  scale_x_continuous(
    breaks = c(56  , 87, 115, 143, 171, 199, 227),
    labels = c(
      "Prime only",
      "4 wk post-boost",
      "8 wk post-boost",
      "12 wk post-boost",
      "16 wk post-boost",
      "20 wk post-boost",
      "24 wk post-boost"
    )
  ) +
  theme(
    axis.text = element_text(size = 12),
    axis.title = element_text(size = 16),
    axis.text.x = element_text(
      angle = 45,
      hjust = 1,
      vjust = 1
    ),
    strip.text = element_text(size = 14),
    legend.text = element_text(size = 14),
    legend.title = element_text(size = 16)
  )

cowplot::plot_grid(
  plotlist = list(raw1, raw2),
  labels = c("A", "B"),
  nrow = 2,
  ncol = 1
)

#ggsave(filename = "figures_tables/ELISA/plot_raw.png", width = 18, height=16)
#ggsave(filename = "figures_tables/ELISA/plot_raw.pdf", width = 18, height=16)

Model fitting

Residual diagnostics when model day as a continuous variable

#-------------------------------------------------------------------------------
# Fit the model
#-------------------------------------------------------------------------------
dtemp <- filter(d, interval != "PBS" & interval != "0") %>%
  mutate(
    dose = factor(dose),
    day = as.numeric(day),
    animal_id = factor(animal_id),
    interval = factor(interval)
  ) %>%
  mutate(
    animal_unique_id = paste0(interval, "_", animal_id),
    animal_unique_id = factor(animal_unique_id)
  )

# fit a gam model
gam_fit <-
  mgcv::gam(
    value ~ interval + dose + s(day, k = 9) + interval:day +
      interval:dose +  dose:day + dose:day:interval +
      s(interval, bs = "re") + s(interval, day, bs = "re") ,
    data = dtemp,
    method = 'REML'
  )


# Residual diagnostics figure
p <- assumption_check(gam_fit)
title <-
  ggdraw() + draw_label(
    "ELISA residual diagnostics when model day as a continuous variable",
    x = 0.03,
    hjust = 0
  )
plot_grid(
  title,
  p,
  ncol = 1,
  rel_heights = c(0.1, 1),
  align = "v"
)

# fit a gam linear model
gam_fit_linear <-
  mgcv::gam(
    value ~ interval + dose + day + interval:day +
      interval:dose +  dose:day + dose:day:interval +
      s(interval, bs = "re") + s(interval, day, bs = "re") ,
    data = dtemp,
    method = 'REML'
  )

# fit a gam model at df = 8
gam_fit_df8 <-
  mgcv::gam(
    value ~ interval + dose + s(day, k = 8) + interval:day +
      interval:dose +  dose:day + dose:day:interval +
      s(interval, bs = "re") + s(interval, day, bs = "re") ,
    data = dtemp,
    method = 'REML'
  )

Diagnostics

The model without smoothing spline on day does not fit as good as the selected model, where the selected model has AIC = 883 and BIC = 1042 and the model without smoothing spline has AIC = 1640 and BIC = 1764.

The selected model has degrees of freedom equals 9 for the smoothing spline, which is the maximum degrees of freedom allowed given the unique covariate combinations. AIC and BIC tests suggested degree of freedom equals 9 has the best fit. The model has df = 8 smoothing spline has AIC = 984 and BIC = 1138.

Therefore, we add smoothing spline with df = 9 to variable day.

# Residual diagnostics figure
p <- assumption_check(gam_fit_linear)
title <-
  ggdraw() + draw_label(
    "ELISA residual diagnostics when model day as a continuous variable",
    x = 0.03,
    hjust = 0
  )
plot_grid(
  title,
  p,
  ncol = 1,
  rel_heights = c(0.1, 1),
  align = "v"
)

Figure 2. S2-P-specific Serum Binding IgG Antibody Titers

#-------------------------------------------------------------------------------
# Figure 2. Spike-specific Serum Binding IgG Antibody Titers
#-------------------------------------------------------------------------------
plot1 <- ggemmeans(gam_fit, terms = c("day", "interval", "dose"))
plot(plot1) +
  scale_color_manual(values = panel0) +
  scale_fill_manual(values = panel0) +
  scale_x_continuous(
    breaks = c(56 , 64, 73, 87, 115, 143, 171, 199, 226),
    labels = c(
      "Prime only",
      "1 wk post-boost",
      "2 wk post-boost",
      "4 wk post-boost",
      "8 wk post-boost",
      "12 wk post-boost",
      "16 wk post-boost",
      "20 wk post-boost",
      "24 wk post-boost"
    )
  ) +
  labs(
    title = expression(Figure2A ~ Predicted ~ values ~ "for" ~ log[10] ~ S2P ~
                         IgG ~ titer),
    y = expression(log[10] ~ S2P ~ IgG ~ titer)
  ) +
  theme(panel.grid.major = element_line(colour = "grey100")) +
  theme(
    axis.text = element_text(size = 12),
    axis.text.x = element_text(
      angle = 45,
      hjust = 1,
      vjust = 1
    ),
    axis.title = element_text(size = 18),
    strip.text = element_text(size = 18),
    legend.text = element_text(size = 16),
    legend.title = element_text(size = 18),
    plot.title = element_text(size = 20),
    legend.position = "bottom"
  ) +
  facet_grid(~ facet)

ggsave(
  filename = "figures_tables/ELISA/ELISA_Titers_F2A.png",
  width = 10,
  height = 8,
  bg = "white"
)
ggsave(
  filename = "figures_tables/ELISA/ELISA_Titers_F2A.pdf",
  width = 12,
  height = 8,
  bg = "white"
)
# plot effect of dose treating dose as numerical variable
plot2 <-
  ggemmeans(
    gam_fit,
    terms = c("day", "dose", "interval"),
    condition = c("day", "interval")
  )
plot(plot2) +
  #scale_color_manual(values=dosepanel) +
  scale_x_continuous(
    breaks = c(56 , 64, 73, 87, 115, 143, 171, 199, 226),
    labels = c(
      "Prime only",
      "1 wk post-boost",
      "2 wk post-boost",
      "4 wk post-boost",
      "8 wk post-boost",
      "12 wk post-boost",
      "16 wk post-boost",
      "20 wk post-boost",
      "24 wk post-boost"
    )
  ) +
  labs(
    x = "",
    title = expression(Figure2B ~ Predicted ~ values ~ "for" ~ log[10] ~ S2P ~
                         IgG ~ titer),
    y = expression(log[10] ~ S2P ~ IgG ~ titer)
  ) +
  theme(panel.grid.major = element_line(colour = "grey100")) +
  theme(axis.text.x = element_text(
    angle = 45,
    hjust = 1,
    vjust = 1
  ),
  legend.position = "bottom") +
  facet_wrap(~ facet)

ggsave(
  filename = "figures_tables/ELISA/ELISA_Titers_F2B.png",
  width = 10,
  height = 8,
  bg = "white"
)
ggsave(
  filename = "figures_tables/ELISA/ELISA_Titers_F2B.pdf",
  width = 12,
  height = 8,
  bg = "white"
)

Residual diagnostics when model day as a factor

# fit the model for mean comparison
dtemp <- filter(d, interval != "PBS" & interval != "0") %>%
  mutate(
    dose = as.factor(dose),
    animal_id = factor(animal_id),
    interval = factor(interval),
    day = factor(day)
  ) %>%
  mutate(animal_id = paste0(interval, "_", animal_id)) %>%
  mutate(animal_id = paste0(animal_id , "_", dose))

lm_fit <-
  lmer(value ~ interval * dose * day + (1 |
                                          animal_id) , data = dtemp)
p2 <- assumption_check(lm_fit)


# Residual diagnostics figure
title <-
  ggdraw() + draw_label(
    "ELISA residual diagnostics when model day as a discrete variable",
    x = 0.03,
    hjust = 0
  )
plot_grid(
  title,
  p2,
  ncol = 1,
  rel_heights = c(0.1, 1),
  align = "v"
)

Supplementary Figure 1. Predicted Fold-change of Spike-specific Serum Binding IgG Antibody Titers

# The results indicate that the peak immune response/ titer levels are achieved at 2 weeks post prime
grid <- ref_grid(lm_fit, cov.keep = c("day", "dose", "interval"))
emm_fit_day <- emmeans(grid, ~ day |  interval + dose)

#-------------------------------------------------------------------------------
# Supplementary Figure 2. Predicted Fold-change of Spike-specific Serum Binding IgG Antibody Titers.
#-------------------------------------------------------------------------------

set.seed(101)
summary(
  contrast(
    emm_fit_day,
    list(
      "2 weeks post-boost \n vs Pre-boost" = c(-1, 0, 1, 0, 0, 0, 0, 0, 0),
      "24 weeks post-boost \n vs Pre-boost" = c(-1, 0, 0, 0, 0, 0, 0, 0, 1)
    )
  ),
  infer = TRUE,
  adjust = "mvt",
  level = .95
)  %>%
  as_tibble() %>%
  mutate(
    .,
    interval = forcats::fct_recode(
      interval,
      "8-week interval" = "8",
      "6-week interval" = "6",
      "4-week interval" = "4",
      "3-week interval" = "3",
      "2-week interval" = "2",
      "1-week interval" = "1"
    )
  ) %>%
  mutate(., interval = factor(.$interval, levels = levels(.$interval)[8:1])) %>%
  mutate(dose = forcats::fct_recode(
    factor(dose),
    "mRNA-1273 1 µg" = "1",
    "mRNA-1273 10 µg" = "10"
  )) %>%
  ggplot(aes(x = contrast , y = estimate, colour = interval)) +
  geom_point(position = position_dodge(width = 0.6), size = 2)  +
  geom_errorbar(
    aes(ymin = lower.CL, ymax = upper.CL),
    width = .1,
    position = position_dodge(width = 0.6),
    size = 1
  ) +
  scale_color_manual(
    values =   c(
      modernared,
      modernablue,
      modernagreen,
      modernapurple,
      modernaorange,
      modernadarkblue2,
      modernadarkgray
    )
  ) +
  theme_bw() +  theme(panel.grid.major = element_line(colour = "grey100")) +
  scale_y_continuous(breaks = log10(c(0.5, 1, 3, 10, 30, 100)), labels = c(0.5, 1, 3, 10, 30, 100)) +
  geom_hline(yintercept = 0,
             color = "darkcyan",
             lty = "dashed") +
  labs(x = "",
       y = expression(Estimated ~ fold ~ change ~ of ~ S2P ~ IgG ~ titer)) +
  facet_grid(~ dose) +
  theme(
    axis.text = element_text(size = 12),
    axis.title = element_text(size = 18),
    strip.text = element_text(size = 18),
    legend.text = element_text(size = 16),
    legend.title = element_text(size = 18),
    plot.title = element_text(size = 20),
    legend.position = "bottom"
  )

ggsave(filename = "figures_tables/ELISA/ELISA_foldchange_SF2.png",
       width = 10,
       height = 8)
ggsave(filename = "figures_tables/ELISA/ELISA_foldchange_SF2.pdf",
       width = 10,
       height = 8)

Supplementary Table 1. Statistical comparisons of S2-P IgG antibody titer in mice between different mRNA-1273 dosing intervals.

#-------------------------------------------------------------------------------
# Supplementary Table 1
#-------------------------------------------------------------------------------
grid <- ref_grid(gam_fit, cov.keep = c("day", "dose", "interval"))
em_fit <- emmeans::emmeans(grid, ~ interval | day + dose)
set.seed(100)
sum_fit <-
  summary(
    pairs(em_fit , reverse = TRUE),
    infer = TRUE,
    adjust = "mvt",
    level = .95
  )


# generate tables
sum_fit %>%
  as_tibble() %>%
  dplyr::select(contrast, estimate, day, dose, p.value) %>%
  mutate(estimate = 10 ^ estimate) %>%
  mutate(contrast = str_replace(contrast, " - ", " / ")) %>%
  rename(
    "Pairwise comparison" = contrast,
    "Fold change" = "estimate",
    "Adjusted P-value" = "p.value"
  ) %>%
  filter(`Adjusted P-value`  < 0.05) %>%
  DT::datatable(caption = "Statistical comparisons of S2-P IgG antibody titer in mice between different mRNA-1273 dosing intervals") %>%
  formatSignif(columns = c('Fold change', 'Adjusted P-value'),
               digits = 7)

Supplementary Figure 4. Statistical comparisons of S2-P IgG Antibody Titers between different mRNA-1273 dosing intervals(Supplementary Table 1 visualization)

sum_fit %>%
  as_tibble() %>%
  mutate(contrast  = str_remove_all(contrast  , "interval")) %>%
  mutate(contrast = str_replace(contrast, " - ", " weeks / ")) %>%
  mutate(contrast = ifelse(
    str_detect(contrast, "1"),
    paste(contrast, "week"),
    paste(contrast, "weeks")
  )) %>%
  mutate(day = as.factor(day)) %>%
  mutate(
    day = forcats::fct_recode(
      factor(day),
      "Prime only" = "56",
      "1 wk post-boost" = "64",
      "2 wk post-boost" = "73",
      "4 wk post-boost" = "87",
      "8 wk post-boost" = "115",
      "12 wk post-boost" = "143",
      "16 wk post-boost" = "171",
      "20 wk post-boost" = "199",
      "24 wk post-boost" = "226"
    )
  ) %>%
  rename(
    "Pairwise comparison" = contrast,
    "Fold change" = "estimate",
    "Adjusted P-value" = "p.value",
    "Day" = "day"
  ) %>%
  mutate(dose = forcats::fct_recode(
    factor(dose),
    "mRNA-1273 1 µg" = "1",
    "mRNA-1273 10 µg" = "10"
  )) %>%
  ggplot(aes(x = `Pairwise comparison` , y = `Fold change`, color = Day)) +
  geom_point(position = position_dodge(width = 2), size = 3)  +
  geom_errorbar(
    aes(ymin = lower.CL, ymax = upper.CL),
    width = .1,
    position = position_dodge(width = 2),
    size = 1
  ) +
  scale_color_manual(
    values = c(
      "#ABCFE5",
      "#93C4DE",
      "#78B5D8",
      "#60A6D1",
      "#4A97C9",
      "#3787C0",
      "#2575B7",
      "#1664AB",
      "#09539D",
      "#084185"
    )
  ) +
  theme_bw() +  theme(panel.grid.major = element_line(colour = "grey100")) +
  scale_y_continuous(breaks = log10(c(0.5, 1, 3, 10, 30, 100)),
                     labels = c(0.5, 1, 3, 10, 30, 100)) +
  geom_hline(yintercept = 0,
             color = "darkcyan",
             lty = "dashed") +
  labs(x = "Prime-boost dosing interval",
       y = expression(Estimated ~ fold ~ change ~ of ~ S2P ~ IgG ~ titer)) +
  facet_grid(`Pairwise comparison` ~ dose,
             scales = "free",
             space = "free") +
  coord_flip() +
  theme(
    axis.text = element_text(size = 18),
    axis.title = element_text(size = 18),
    strip.text = element_text(size = 18),
    legend.text = element_text(size = 16,
                               margin = margin(t = 10)),
    legend.title = element_text(size = 18),
    plot.title = element_text(size = 20),
    strip.text.y = element_blank(),
    legend.position = "right"
  )

ggsave(filename = "figures_tables/ELISA/ELISA_emmeans.png",
       width = 16,
       height = 18)
ggsave(filename = "figures_tables/ELISA/ELISA_emmeans.pdf",
       width = 16,
       height = 18)

Session information

sessionInfo()
## R version 3.6.3 (2020-02-29)
## Platform: x86_64-apple-darwin15.6.0 (64-bit)
## Running under: macOS Catalina 10.15.7
## 
## Matrix products: default
## BLAS:   /Library/Frameworks/R.framework/Versions/3.6/Resources/lib/libRblas.0.dylib
## LAPACK: /Library/Frameworks/R.framework/Versions/3.6/Resources/lib/libRlapack.dylib
## 
## locale:
## [1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8
## 
## attached base packages:
## [1] stats     graphics  grDevices utils     datasets  methods   base     
## 
## other attached packages:
##  [1] cowplot_1.1.1     DT_0.24           pzfx_0.3.0        ggpubr_0.4.0     
##  [5] formattable_0.2.1 lme4_1.1-28       Matrix_1.2-18     ggeffects_1.1.3  
##  [9] emmeans_1.7.5     mgcv_1.8-31       nlme_3.1-144      readxl_1.4.0     
## [13] forcats_0.5.1     stringr_1.4.0     dplyr_1.0.9       purrr_0.3.4      
## [17] readr_2.1.2       tidyr_1.2.0       tibble_3.1.8      ggplot2_3.3.6    
## [21] tidyverse_1.3.2  
## 
## loaded via a namespace (and not attached):
##  [1] pbkrtest_0.5.1      fs_1.5.2            lubridate_1.8.0    
##  [4] insight_0.18.2      httr_1.4.3          tools_3.6.3        
##  [7] backports_1.4.1     bslib_0.4.0         sjlabelled_1.2.0   
## [10] utf8_1.2.2          R6_2.5.1            DBI_1.1.3          
## [13] colorspace_2.0-3    withr_2.5.0         tidyselect_1.1.2   
## [16] compiler_3.6.3      textshaping_0.3.6   cli_3.3.0          
## [19] rvest_1.0.2         xml2_1.3.3          labeling_0.4.2     
## [22] sass_0.4.2          scales_1.2.0        mvtnorm_1.1-3      
## [25] systemfonts_1.0.2   digest_0.6.29       minqa_1.2.4        
## [28] rmarkdown_2.14      pkgconfig_2.0.3     htmltools_0.5.3    
## [31] highr_0.9           dbplyr_2.2.1        fastmap_1.1.0      
## [34] htmlwidgets_1.5.4   rlang_1.0.4         rstudioapi_0.13    
## [37] farver_2.1.1        jquerylib_0.1.4     generics_0.1.3     
## [40] jsonlite_1.8.0      crosstalk_1.2.0     car_3.1-0          
## [43] googlesheets4_1.0.1 magrittr_2.0.3      Rcpp_1.0.9         
## [46] munsell_0.5.0       fansi_1.0.3         abind_1.4-5        
## [49] lifecycle_1.0.1     stringi_1.7.8       yaml_2.3.5         
## [52] carData_3.0-5       MASS_7.3-51.5       grid_3.6.3         
## [55] parallel_3.6.3      crayon_1.5.1        lattice_0.20-38    
## [58] haven_2.5.0         splines_3.6.3       hms_1.1.1          
## [61] knitr_1.39          pillar_1.8.0        boot_1.3-24        
## [64] estimability_1.4.1  ggsignif_0.6.3      reprex_2.0.1       
## [67] glue_1.6.2          evaluate_0.16       modelr_0.1.8       
## [70] vctrs_0.4.1         nloptr_1.2.2.3      tzdb_0.3.0         
## [73] cellranger_1.1.0    gtable_0.3.0        assertthat_0.2.1   
## [76] cachem_1.0.6        xfun_0.32           xtable_1.8-4       
## [79] broom_1.0.0         coda_0.19-4         rstatix_0.7.0      
## [82] ragg_1.1.3          googledrive_2.0.0   gargle_1.2.0       
## [85] ellipsis_0.3.2